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A theory is presented which allows us to accurately calculate the density profile of monovalent 
and multivalent counterions in suspensions of polarizable colloids or nano-particles. In the case 
of monovalent ions, we derive a weak-coupling theory that explicitly accounts for the ion-image 
interaction, leading to a modified Poisson-Boltzmann equation. For suspensions with multivalent 
counterions, a strong-coupling theory is used to calculate the density profile near the colloidal 
surface and a Poisson-Boltzmann equation with a renormalized boundary condition to account for 
the counterion distribution in the far-field. All the results are compared with the Monte Carlo 
simulations, showing an excellent agreement between the theory and the simulations. 
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Colloidal suspensions are of great practical interest for 
biology, chemistry, and physics. The subject has a long 
history going back more than a hundred years. In spite 
of the intense effort, many interesting phenomena which 
are found in colloidal science have not been fully eluci- 
dated. For example, it is well known that the stability 
of a hydrophobic colloidal suspension depends specifically 
on the electrolyte present in suspension. Addition of mul- 
tivalent counterions results in a rapid precipitation of col- 
loidal particles. What is more surprising is that even for 
monovalent counterions, stability of colloidal suspensions 
depend strongly on the precise nature of the counteri- 
ons. Thus, as one goes along the halogen series, the crit- 
ical coagulation concentrations of positive colloidal par- 
ticles can decrease by as much as an order of magnitude, 
when anion is changed from fluoride to iodide l]. An- 
other interesting phenomenon found in suspensions with 
multivalent ions is the reversal of electrophoretic mobil- 
ity @, H|, or equivalently charge reversal 0-0] • Under 
some conditions, it is also possible to observe like-charge 
attraction between the colloidal particles of the same sign 
of charge [§-12|. Many of these interesting phenomena 



are the consequence of strong electrostatic correlations 
between the counterions. The role of electrostatic cor- 
relations has been studied using simple models of col- 
loidal suspensions which neglect particle polarizability. 
The standard Poisson-Boltzmann equation (PB) — used 
extensively in colloidal science — fails to account for the 
induced charge at the particle-solvent interface, predict- 
ing that the counterion density should remain unaffected 
by the colloidal polarizability. In this Letter, we will show 
that the induced colloidal charge significantly modifies 
the ionic density distribution even for monovalent coun- 
terions. The theory developed in this Letter allows us 



to accurately predict the counterion density distribution 
both in the weak (monovalent counterions) and strong 
(multivalent counterions) coupling limits. 

We will use the primitive model of colloidal suspen- 
sion in which colloidal particles are represented by hard 
spheres of radius a and dielectric constant e c with the 
charge —Zq distributed uniformly over the surface. Wa- 
ter is modeled as a uniform dielectric of permittivity e w . 
The system is at room temperature, so that the Bjerrum 
length, defined as \b = q 2 / e w ksT , is 7.2A. To account 
for the finite concentration of colloidal particles, we will 
use a spherical Wigner-Seitz (WS) cell of radius R. The 
cell also contains N — Z/a a-valent counterions, each of 
radius r c . For most colloidal suspensions of practical in- 
terest e w /e c 3> 1. In this limit it is possible to show that 
the exact Green function for the interaction between two 
counterions is 11311. 
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The terms in Eq. [T] are respectively the electrostatic po- 
tential produced by an ion located at r', image charge 

2 

located at the inversion point pj-r' inside the colloid, 
and the counter-image charge spread uniformly from the 
inversion point up to the center of the colloidal parti- 
cle dlGlj, 
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Eq. [T] can be used to obtain the counterion density dis- 
tributions using Monte Carlo (MC) as was shown in Ref. 
[13I ] , and is a much faster alternative to simulations based 
on expansion in Legendre polynomials fl6j ] . 
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To calculate the density profiles theoretically one of- 
ten relies on the mean-field PB equation. In the case of 
non-polarizablc colloidal particles, it is well known that 
for monovalent counterions the PB theory is very accu- 
rate 17j . On the other hand, for multivalent counterions 
there are str ong deviations from the prediction of PB 
equation 17H22|. These deviations are a consequence of 



strong electrostatic correlation between the multivalent 
counterions close to the colloidal surface. For polarizable 
particles, however, the usual PB equation fails even for 
monovalent counterions. In this case, the failure of the 
PB equation is a consequence of the counterion-image 
interaction — i.e. transverse correlations — which are 
left out of the mean-field PB equation. The deviations 
from PB theory augment with increasing counterion va- 
lence. One of the objectives of the present work is to 
construct a theory that can properly take into account 
both the ion-ion and the ion-image interactions for po- 
larizable colloidal particles. 

The electrostatic potential inside the WS cell satisfies 
the exact Poisson equation 
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where 4>{r) is the mean electrostatic potential at distance 
r from the colloidal center and p(r) is the mean counte- 
rion density. We next suppose that the counterions are 
distributed in accordance with the Boltzmann distribu- 
tion 
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where (3 = l/fegT and z is the distance from the col- 
loidal surface, z = r — a. In writing Eq. |4l we have 
implicitly assumed that the main contribution to the po- 
tential of mean-force, besides the electrostatic potential 
4>(r), comes from the ion-image interaction, W(k, z). Al- 
though this is true for the monovalent counterions - 
weak-coupling limit — this assumption breaks-down for 
the multivalent ions in the strong-coupling limit. 

Unfortunately there is no exact way of calculating the 
charge-image interaction, so that approximations must 
be invoked. We proceed as follows. Consider a neutral 
one component plasma (OCP) confined to a half-space by 
a hard dielectric wall of e c m 0. For z < r c the hard core 
repulsion requires that the electrostatic potential satis- 
fies the Laplace equation V 2 ^(r) = 0, where we have 
assumed that the neutralizing background of the OCP 
also starts at z — r c . For z > r c the electrostatic poten- 
tial inside the OCP satisfies the linearized PB equation 

V 2 0(r) = K 2 4>{r), where k = \J iva y BZ and V is the vol- 
ume accessible to the counterions. These two equations 
can be solved explicitly to calculate the work that must 
be done to bring an ion of the OCP from infinity to the 



distance z from the wall 1 23.1. We obtain 



(5) 



where 
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and s — \/k 2 + k 2 . Eq. [5] accounts for two funda- 
mental contributions: the ion-image interaction which 
is screened with the Debye length 1/k, and for the loss 
of the electrostatic solvation free energy experienced by 
the ion near the interface. Since the wall has the dielec- 
tric constant much smaller than that of water, the ion- 
image interaction is strongly repulsive. However, note 
that even if there would not be any dielectric discon- 
tinuity but simply a hard wall, there still would be an 
electrostatic energy penalty for bringing an ion from the 
bulk of electrolyte — where its self-energy is efficiently 
screened by the background and by the other ions - 
to the interface, where the broken symmetry prevents 
the efficient screening of its electric field, resulting in a 
higher electrostatic self-energy. Therefore, besides the 
ion-image repulsion, at finite ionic concentrations there 
is an additional repulsive interaction arising from the bro- 
ken translational symmetry near the wall- water interface, 
which is also taken into account in Eq. [5j 
It is instructive to rewrite Eq. [5] as 
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The first term in the square brackets, A sc (k, z) = 
/3W ( K ,r c )r c /z - /3W (0, r c )r c /z with /3W (0,r c ) = 
a 2 As/4r c , is the repulsive energy arising from the bro- 
ken translational symmetry, while the second term is 
due to the ion-image interaction. For colloids and nano- 
particles, we expect that A sc (k, z) depends only weakly 
on the curvature and can be approximated by that of a 
planar wall. On the other hand, the ion-image interaction 
energy is strongly dependent on the radius of curvature of 
the dielectric interface and must be corrected when treat- 
ing nano-particles or colloids. Using the Green function 
Eq. [T] the curvature correction can be easily calculated. 
We find 



0W(k,z) 
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These expressions are used in Eqs.|3]and|2]to account for 
the colloidal polarizability. The modified PB equation is 
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then solved numerically to calculate the counterion den- 
sity profiles. In Fig. [T] the theory is compared with the 
MC simulations for monovalent ions, showing an excel- 
lent agreement. On the other hand, in the Supplemen- 
tary Information, we show that the usual PB equation 
deviates significantly from the MC data near the colloidal 
surface. 




■ Z= 1500\ 

100 200 300 400 150 300 450 600 

r[A] r[A] 



FIG. 1. Density profiles of monovalent counterions. Symbols 
represent simulation data and lines represent the modified 
PB theory. The parameters of the simulations are a = 1, 
r c = 2A; down triangles (a = 150A and R = 400A); up 
triangles (a = 100A and R = 400A); diamonds (a = 150A 
and R = 600A). The insets show the profiles close to the 
colloidal surface. 

Although we find a perfect agreement between the 
theory and the simulations for monovalent counterions, 
strong deviations are observed when multivalent ions are 



present in suspension 241] . This is similar to what has 



been previously found with non-polarizable colloidal par- 
ticles [2(| . Near the colloidal surface, ionic concentration 
is very large and the electrostatic interaction between the 
multivalent counterions is very strong. This results in a 
formation of a strongly correlated quasi-two-dimensional 
OCP on the surface of colloidal particle. The strength of 
electrostatic correlations can be quantified by the plasma 



parameter T 



a 3/2 \ B \/Z 
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which is the ratio of the char- 



acteristic electrostatic to the thermal energy of condensed 
counterions For r ^> 1 the counterions exhibit local 
hexagonal order. 

Consider a counterion of the condensed layer. If the 
average separation between the condensed counterions is 
much smaller than the radius of the nano-particle, the ef- 
fects of the curvature will be screened. Thus, in the limit 



r > 1, a > r c , and 



<C 1, we can neglect the cur- 



(a+r c )r 

vature of the colloidal surface and treat it as a hard wall 
separating the low dielectric e c « half-space from the 
high-dielectric region occupied by the counterions, Fig. [2] 
As a consequence of the hexagonal symmetry, the elec- 
tric field produced on the counterion by other condensed 
counterions vanishes. Therefore, in the limit T — ¥ oo, the 
ion interacts only with the electric field of the charged 
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FIG. 2. Hexagon of images at the surface. In (A) the side 
view. In (B) the self-image and the nearest neighbors. 



wall and the field produced by the images. The average 
separation h between the condensed counterions of the 

hexagonal lattice is ft, = \J 8ir(a + r c ) 2 a/ (Zy/3), consid- 
ering that all Z/a are in the area 47r(a + r c ) 2 . In the 
strong coupling limit, the counterions density distribu- 
tion then take a particularly simple form 



p sc {z) = A e -2(*-rc)/; 9 e-/W„(z) 



(10) 



where l gc = 1S ^ ne Gouy- Chapman length and 

Wh{z) is the electrostatic potential of interaction between 
the ion and its image and the images of its first nearest 
neighbors, 



PWk(z) = 
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The first term in the exponential of Eq. [TU] is due to the 
ion interacting with the field of the charged hard wall. 
Note that because e c ~ 0, the displacement field inside 
the low-dielectric half-space vanishes, which leads to the 
factor of 2 in the exponent of Eq. [101 that is absent in the 
strong-coupling theories of non-polarizable particles [20j j . 
The normalization factor A is obtained using the charge 
neutrality condition, 
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In Fig. |3] we show the density distributions for different 
colloidal systems containing either divalent or trivalent 
counterions, a — 2 and a — 3, respectively. A very good 
agreement is observed between the theory and the simu- 
lations. As expected, the agreement improves for larger 
r and when a — > oo (see Supplementary Information). 

Far from the colloidal surface the concentration of 
counterions drops rapidly. In this far-field region the 
counterion correlations can be neglected and the mean- 
field Poisson-Boltzmann theory once again becomes ap- 
plicable. To connect the strongly-correlated region near 
the colloidal surface with the weakly-correlated far-field, 
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we use the theory developed in Ref . [20j ■ The uniformity 
of the chemical potential throughout the system requires 
that the countcrion concentration at the boundary of the 
WS cell be related with the coarse-grained density in the 
strongly-correlated region described by Eq. [TU] resulting 
in equation 



Pbc{R) = Pc g e 



Paq[<i>(a+r a )-<f>{R)}+pn c 



(13) 



where [3p c = -1.65r + 2.6LF 1 / 4 - 0.26 lnT - 1.95 is the 
chemical potential of the strongly-correlated 2d OCP 25] 
and p cg is the coarse-grained density of the 2d OCP, 



fr° +S " dz p sc (z) 



(14) 



where l sc — 3.6 l gc [20(. Solving the usual PB equation 
while enforcing the condition Eq. [13] at the cell bound- 
ary, we obtain both the effective colloidal charge and the 
countcrion density profiles in the weakly-correlated re- 
gion. In Fig. HI we show the comparison between MC 
simulations and the predictions of the present theory in 
the far-field region. Once again an excellent agreement 
is found between the theory and the simulations. 

We have derived a new theory which allows us to quan- 
titatively calculate the density distributions of monova- 
lent and multivalent counterions in suspensions of polar- 
izable colloidal particles. In the case of the monovalent 
ions, we have derived a modified PB equation which al- 
lows us to calculate the complete countcrion density pro- 
file. For multivalent ions, we used the strong-coupling 
theory to obtain the density profiles near the colloidal 
surface and the usual PB equation with the renormalized 
boundary condition to calculate the density distribution 
in the far field. All the results are in excellent agree- 
ment with the MC simulations. Finally, we mention that 
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FIG. 3. Density profiles near the colloidal surface. Symbols 
are the simulation data and lines are the predictions of the 
present theory. Circles, squares and diamonds are simulation 
data for a = 3, a = 150A and R — 600A. Right, left triangles 
and stars are simulation data for a = 2, a — lOOA and R = 
400 A. The ionic radius is r c = 2 A. 
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FIG. 4. Density profiles in the far-field. Symbols represent 
the simulation data and lines are the predictions of the present 
theory. The parameters are the same as in Fig. [3] 



the same theory can be easily extended to treat metal 
nano-particles. 
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